library(lfe)
library(lme4)
library(foreign)
library(DataCombine)
library(plyr)
library(survival)
library(Zelig)
library(doBy)
library(interplot)
library(mgcv)
library(stargazer)
library(maps)
library(maptools)
library(ggmap)
library(mapdata)
library(RgoogleMaps)
#Read data
setwd("C:/Projects/Replication files/Data/")
fs.dat <- read.dta("fs.dat6.dta")
###Correlations between maize and urbanization and maize and nighttime light
#Urbanization
maize_urb <- na.omit(fs.dat[,c("lagurban", "maize_area")])
cor(maize_urb$lagurban, maize_urb$maize_area)
#Nighttimelight
maize_nl <- na.omit(fs.dat[,c("nlights_calib_mean", "maize_area")])
cor(maize_nl$nlights_calib_mean, maize_nl$maize_area)

#Number of countries with maize that experienced conflict
bb <- summaryBy(maize_area + civconf ~ ccode, FUN=c("sum"), keep.names = T, data=fs.dat, na.rm=T)
bb$maize_ind <- ifelse(bb$maize_area>0,1,0)
bb$civconf <- ifelse(bb$civconf>0,1,0)
table(bb$civconf, bb$maize_ind)
(58/(58+2))*100
#Cells and average cell years with conflict
fs.dat$one <- 1
col <- summaryBy(one ~ gid + ccode, data=fs.dat, keep.names = TRUE)
col$one <- 1
col.s <- summaryBy(one~ccode, data=col, FUN=c("sum"), keep.names = TRUE)
mean(col.s$one)
median(col.s$one)

###################
###################
###Plot the maps###
###################
###################

##########################
###Rebellion (Figure 1)###
##########################
dat.2 <- read.dta("reb.map.dat2.dta")
###Create data for plotting
#Subset out
dat.2.s <- na.omit(dat.2[,c("latitude", "longitude", "civconf")])
values2 <- summaryBy(civconf ~ latitude + longitude, data=dat.2.s, FUN=c("sum"), keep.names = TRUE, na.rm=TRUE)
#values2 <- na.omit(dat.2[,c("latitude", "longitude", "civconf")])
values$civconf <- values$civconf/(max(dat.2$year)-min(dat.2$year)+1)
values2$conf_norm <- values2$civconf/max(values2$civconf)
#Create a color spectrum for the map
col.s2 <- rgb(1, 0, 0, alpha=1*values2$conf_norm, maxColorValue = 1)
#Create color spectrum for the legend
c02 <- rgb(1, 1, 1, maxColorValue = 1)
c12 <- rgb(1, 0, 0, alpha=0.33, maxColorValue = 1)
c22 <- rgb(1, 0, 0, alpha=0.66, maxColorValue = 1)
c32 <- rgb(1, 0, 0, alpha=1, maxColorValue = 1)
##Vector of categories
colvec.at2 <- rbind(c02, c12, c22, c32)
#Create categories for the legend
leg.at.txt2 <- c("No rebellion", "Low intensity", "Medium instensity", "High intensity")
###Plot the map
jpeg("confmap.jpeg", width = 7, height = 6, units = 'in', res = 500)
map("worldHires", col="white", fill=TRUE,xlim=c(-180, 180), ylim=c(-58, 90))
points(values2$longitude, values2$latitude, pch=16, col=col.s2, cex=0.15)
legend("bottomleft",title = "Average Rebellion Intensity by 0.5 Degree Grids",
       legend = leg.at.txt2,
       fill=colvec.at2,
       cex = 0.56,
       bty = "n")
dev.off()

###########################
###Maize area (Figure 2)###
###########################
dat.1 <- read.dta("maize.map.dat2.dta")
dat.1.s <- na.omit(dat.1[,c("latitude", "longitude", "maize_area")])
values <- summaryBy(maize_area ~ latitude + longitude, FUN=c("sum"), data=dat.1.s, keep.names = TRUE, na.rm=TRUE)
values$maize_area <- values$maize_area/(max(dat.1$year)-min(dat.1$year)+1)
###Create data for plotting
#Subset out
#values$lagmaize_area <- ifelse(values$lagmaize_area > 0, values$lagmaize_area+0.0001, 0)
values$maize_norm <- values$maize_area/max(values$maize_area)
#Create a color spectrum for the map
col.s <- rgb(0, 0, 1, alpha=1*values$maize_norm, maxColorValue = 1)
#Create color spectrum for the legend
c0 <- rgb(1, 1, 1, maxColorValue = 1)
c1 <- rgb(0, 0, 1, alpha=0.33, maxColorValue = 1)
c2 <- rgb(0, 0, 1, alpha=0.66, maxColorValue = 1)
c3 <- rgb(0, 0, 1, alpha=1, maxColorValue = 1)
##Vector of categories
colvec.at <- rbind(c0, c1, c2, c3)
#Create categories for the legend
leg.at.txt <- c("No maize", "Low productivity", "Medium productivity", "High productivity")
###Plot the map
jpeg("maizemap.jpeg", width = 7, height = 6, units = 'in', res = 500)
map("worldHires", col="white", fill=TRUE,xlim=c(-180, 180), ylim=c(-58, 90))
points(values$longitude, values$latitude, pch=15, col=col.s, cex=0.3)
legend("bottomleft",title = "Average Maize Area by 0.5 Degree Grids",
       legend = leg.at.txt,
       fill=colvec.at,
       cex = 0.56,
       bty = "n")
dev.off()

#####################################################
#####################################################
###Main Analysis (Tables 1 and 2, Figures 3 and 4)###
#####################################################
#####################################################

####################
###Baseline model###
####################
###Subset NA omit data
ss.dat.b <- na.omit(fs.dat[,c("civconf", "pop", "maize_area", "nlights_calib_mean", "lagcivconflagtemp", "gid", "year")])
ss.dat.b$logpop <- log(ss.dat.b$pop+1)
#Create interation
ss.dat.b$light_main <- ss.dat.b$nlights_calib_mean*ss.dat.b$maize_area

###Run regular GLM
glm.mod.s <- glm(civconf ~ maize_area+nlights_calib_mean+light_main+
                   logpop+lagcivconflagtemp+
                   factor(year)
                 ,data=ss.dat.b, family=binomial(link = "logit"))
summary(glm.mod.s)

###Run GAM
gam.mod.s2 <- bam(civconf ~ s(maize_area, bs='tp', sp=0.2) +
                    s(nlights_calib_mean, bs='tp', sp=0.2) +
                    ti(nlights_calib_mean, maize_area, bs=c('tp','tp'), sp=c(0.2,0.2)) +
                    s(pop, bs='tp', sp=0.2)+
                    lagcivconflagtemp+
                    factor(year)
                  ,data=ss.dat.b, family=binomial(link = "logit"))
summary(gam.mod.s2)
gam.mod.s2$aic
##Plot diagnostics (Figure A2)
jpeg("diagbase2.jpeg", width = 6, height = 6, units = 'in', res = 500)
par(mfrow=c(2,2))
plot(gam.mod.s2,se = TRUE)
dev.off()

################
###Full model###
################
###Subset NA omit data
ss.dat.f <- na.omit(fs.dat[,c("civconf", "drug_y", "splcivconf", "nwstate", "temp", "spi6", "logttime", "pop", "polity2", "maize_area", "nlights_calib_mean", "lagcivconflagtemp", "gid", "year")])
ss.dat.f$logpop <- log(ss.dat.f$pop+1)
ss.dat.f$ttime <- expm1(ss.dat.f$logttime)
#Create interaction term
ss.dat.f$light_main <- ss.dat.f$nlights_calib_mean*ss.dat.f$maize_area

###Run regular GLM
glm.mod.f <- glm(civconf ~ maize_area+nlights_calib_mean+light_main+
                   logpop+lagcivconflagtemp+polity2+
                   logttime+temp+spi6+drug_y+splcivconf+nwstate+
                   factor(year)
                 ,data=ss.dat.f, family=binomial(link = "logit"))
summary(glm.mod.f)

###Plot the interaction (Figure 3)
#Model
glm.mod.fi <- glm(civconf ~ maize_area+nlights_calib_mean+
                    maize_area*nlights_calib_mean+
                    logpop+polity2+
                    logttime+temp+spi6+drug_y+splcivconf+nwstate+
                    lagcivconflagtemp+
                    factor(year)
                  ,data=ss.dat.f, family=binomial(link = "logit"))
##Plot Figure 3
jpeg("intglmfull.jpeg", width = 6, height = 6, units = 'in', res = 500)
interplot(m = glm.mod.fi, var1 = "nlights_calib_mean", var2 = "maize_area")+
  # Add labels for X and Y axes
  xlab("Maize") +
  ylab("Estimated Coefficient for\nNighttime Light") +
  # Change the background
  theme_bw() +
  # Add the title
  #ggtitle("Estimated Effect of Maize \non Nighttime Light's Propensity of Rebellion \n(Full Model)") +
  theme(plot.title = element_text(face="bold")) +
  # Add a horizontal line at y = 0
  geom_hline(yintercept = 0, linetype = "dashed")
dev.off()

###Run GAM
gam.mod.f2 <- bam(civconf ~ s(maize_area, bs='tp', sp=0.2) +
                    s(nlights_calib_mean, bs='tp', sp=0.2) +
                    ti(nlights_calib_mean, maize_area, bs=c('tp','tp'), sp=c(0.2,0.2)) +
                    s(pop, bs='tp', sp=0.2)+
                    lagcivconflagtemp+s(ttime, bs='tp', sp=0.2) +
                    polity2 + temp + spi6 + drug_y + splcivconf + nwstate+
                    factor(year)
                  ,data=ss.dat.f, family=binomial(link = "logit"))
summary(gam.mod.f2)
gam.mod.f2$aic

##Plot diagnostics (Figure A3)
jpeg("diagfull2.jpeg", width = 6, height = 6, units = 'in', res = 500)
par(mfrow=c(3,2))
plot(gam.mod.f2,se = TRUE)
dev.off()

##Plot Figure 4
jpeg("intgamfull2.jpeg", width = 6, height = 6, units = 'in', res = 500)
crop.f <- with(ss.dat.f, seq(min(maize_area), max(maize_area), len=100))
nl.f <- with(ss.dat.f, seq(min(nlights_calib_mean), max(nlights_calib_mean), len=100))
new.data.f <- expand.grid(nlights_calib_mean=nl.f, maize_area=crop.f, pop=mean(ss.dat.f$pop), lagcivconflagtemp=0, ttime=mean(ss.dat.f$ttime), polity2=median(ss.dat.f$polity2), temp=mean(ss.dat.f$temp), spi6=median(ss.dat.f$spi6), drug_y=0, splcivconf=0, nwstate=0, year=2000)
pred.at.f <- matrix(predict.bam(gam.mod.f2, new.data.f, type="response"), 100, 100, 100)
persp(crop.f, nl.f, pred.at.f, theta=45, phi=30, ticktype = "detailed",
      xlab="Maize Area", ylab= "Nighttime Light", zlab="Conflict",
      shade =0.5, nticks = 4)
dev.off()


###Export GLMs to LaTex
stargazer(glm.mod.s, glm.mod.f)


#####################
#####################
###Online Appendix###
#####################
#####################

######################################
###Robustness GLM Models (Table A3)###
######################################
#Create interation
fs.dat$light_main <- fs.dat$nlights_calib_mean*fs.dat$maize_area

###Model A1: Added geospatial controls
glm.mod.f.geo <- glm(civconf ~ maize_area+nlights_calib_mean+light_main+
                       logpop+lagcivconflagtemp+polity2+
                       logttime+temp+spi6+drug_y+splcivconf+nwstate+
                       mnt1 + logbdist1 + logcellarea +
                       factor(year)
                     ,data=fs.dat, family=binomial(link = "logit"))
summary(glm.mod.f.geo)
AIC(glm.mod.f.geo)

###Model A2: Geospatial controls, instability, fractionalization, and natural resources
glm.mod.f.cont <- glm(civconf ~ maize_area+nlights_calib_mean+light_main+
                        logpop+lagcivconflagtemp+polity2+
                        logttime+temp+spi6+drug_y+splcivconf+nwstate+
                        mnt1 + logbdist1 + logcellarea +
                        instability + ethfrac + relfrac + gem_y +
                        logross_oil_prod + logross_gas_prod +
                        factor(year)
                      ,data=fs.dat, family=binomial(link = "logit"))
summary(glm.mod.f.cont)
AIC(glm.mod.f.cont)

###Model A3: CSEs
glm.mod.f.cse <- glm(civconf ~ maize_area+nlights_calib_mean+light_main+
                       logpop+lagcivconflagtemp+polity2+
                       logttime+temp+spi6+drug_y+splcivconf+nwstate+
                       factor(year) + cluster(gid)
                     ,data=fs.dat, family=binomial(link = "logit"))
summary(glm.mod.f.cse)
AIC(glm.mod.f.cse)

###Model A4: CFEs
glm.mod.f.cfe <- glm(civconf ~ maize_area+nlights_calib_mean+light_main+
                       logpop+lagcivconflagtemp+polity2+
                       logttime+temp+spi6+drug_y+splcivconf+nwstate+
                       factor(year) + factor(ccode)
                     ,data=fs.dat, family=binomial(link = "logit"))
summary(glm.mod.f.cfe)
AIC(glm.mod.f.cfe)

###Model A5: Country random effects logit
glm.mod.f.re <- glmer(civconf ~ maize_area+nlights_calib_mean+light_main+
                        logpop+lagcivconflagtemp+polity2+
                        logttime+temp+spi6+drug_y+splcivconf+nwstate+
                        (1|year) + (1|ccode),
                      data=fs.dat, family=binomial(link="logit"))
summary(glm.mod.f.re)
#AIC(glm.mod.f.re)

####Model A6: Rare events
#Model
glm.mod.f.rare <- zelig(civconf ~ maize_area+nlights_calib_mean+light_main+
                          logpop+lagcivconflagtemp+polity2+
                          logttime+temp+spi6+drug_y+splcivconf+nwstate+
                          factor(year),
                        data=ss.dat.f, model = "relogit")
summary(glm.mod.f.rare)

###############################
###GAM Robustness (Table A4)###
###############################
###Run GAM baseline
gam.mod.s3 <- bam(civconf ~ s(maize_area,   bs='ps', sp=0.2) +
                    s(nlights_calib_mean,   bs='ps', sp=0.2) +
                    ti(nlights_calib_mean, maize_area, bs=c('ps','ps'), sp=c(0.2,0.2)) +   
                    s(pop,  bs='ps', sp=0.2)+ lagcivconflagtemp +
                    factor(year)
                  ,data=ss.dat.b, family=binomial(link = "logit"))
summary(gam.mod.s3)
gam.mod.s3$aic
##Plot Diagnostics (Figure A4)
jpeg("diagbase3.jpeg", width = 6, height = 6, units = 'in', res = 500)
par(mfrow=c(2,2))
plot(gam.mod.s3,se = TRUE)
dev.off()


###Run GAM full
gam.mod.f3 <- bam(civconf ~ s(maize_area, bs='ps', sp=0.2) +
                    s(nlights_calib_mean, bs='ps', sp=0.2) +
                    ti(nlights_calib_mean, maize_area, bs=c('ps','ps'), sp=c(0.2,0.2)) +
                    s(pop, bs='ps', sp=0.2)+
                    lagcivconflagtemp+s(ttime, bs='ps', sp=0.2) +
                    polity2 + temp + spi6 + drug_y + splcivconf + nwstate+
                    factor(year)
                  ,data=ss.dat.f, family=binomial(link = "logit"))
summary(gam.mod.f3)
gam.mod.f3$aic
##Plot Diagnostics (Figure A5)
jpeg("diagfull3.jpeg", width = 6, height = 6, units = 'in', res = 500)
par(mfrow=c(3,2))
plot(gam.mod.f3,se = TRUE)
dev.off()


#####################################################
###Maize and Cropland Correlation Plot (Figure A1)###
#####################################################

#Import data
b <- read.dta("maize.crop.dta")
#Correlations
cor(b$CropAreaFraction, b$lagmaize_area)
#Trendline
lm.pred <- lm(CropAreaFraction ~ lagmaize_area, data=b)
#Confidence intervals
pp.pred <- predict(lm.pred)
V <- vcov(lm.pred)
X <- model.matrix(~ lagmaize_area, data=b)
se.fit <- sqrt(diag(X %*% V %*% t(X)))
predframe <- with(b, data.frame(lagmaize_area,
                                CropAreaFraction=pp.pred,lwr=pp.pred-1.96*se.fit,upr=pp.pred+1.96*se.fit))

p <- ggplot(b, aes(x = lagmaize_area, y = CropAreaFraction))
#Plot Figure A1
pdf("maizecropcor.pdf")
p +  labs(x="Average Maize % Cell Area (1993-2009)", y="All Cropland (2000)") +
  geom_point()+ ylim(0,1) + xlim(0, 0.6) +
  annotate("text", x = 0.55, y = 0.95, label = "Correlation = 0.41") +
  geom_line(data=predframe, size=0.5, color="green") +
  geom_ribbon(data=predframe,aes(ymin=lwr,ymax=upr),alpha=0.3)
dev.off()

####################################
###Country Plots (Figures A8-A13)###
####################################
###Subset na.omit data
state.dat <- na.omit(fs.dat[,c("ccode", "year", "gid", "maize_area", "nlights_calib_mean", "civconf", "latitude", "longitude")])

###Burundi
##Correlations
#Subset
bur.dat <- subset(state.dat, ccode==516)

##Maps (Figure A8)
#Collapse data
values.b1 <- summaryBy(maize_area ~ latitude + longitude, keep.names = TRUE, FUN=c("mean"), data=bur.dat, na.rm=TRUE)
values.b2 <- summaryBy(nlights_calib_mean ~ latitude + longitude, keep.names = TRUE, FUN=c("mean"), data=bur.dat, na.rm=TRUE)
values.b3 <- summaryBy(civconf ~ latitude + longitude, keep.names = TRUE, FUN=c("mean"), data=bur.dat, na.rm=TRUE)
#Color patterns
col.s1 <- rgb(0,1,0, alpha=1*(values.b1$maize_area/max(values.b1$maize_area)), maxColorValue=1)
col.s2 <- rgb(0,0,1, alpha=1*(values.b2$nlights_calib_mean/max(values.b2$nlights_calib_mean)), maxColorValue=1)
col.s3 <- rgb(1,0,0, alpha=1*(values.b3$civconf/max(values.b3$civconf)), maxColorValue=1)
##Plot
pdf("burudnimap.pdf")
#Split to three columns
par(mfrow=c(1,3))
#Maize
map("worldHires", "Burundi", col="white", fill=TRUE)
points(values.b1$longitude, values.b1$latitude, pch=15, col=col.s1, cex=10)
#Nighttime light
map("worldHires", "Burundi", col="white", fill=TRUE)
points(values.b2$longitude, values.b2$latitude, pch=15, col=col.s2, cex=10)
#Rebellion
map("worldHires", "Burundi", col="white", fill=TRUE)
points(values.b3$longitude, values.b3$latitude, pch=15, col=col.s3, cex=10)
dev.off()

###Correlations (Figure A9)
cor(bur.dat$civconf, bur.dat$maize_area)
cor(bur.dat$civconf, bur.dat$nlights_calib_mean)
pdf("burmaize.pdf")
p <- ggplot(bur.dat, aes(x = maize_area, y = civconf))
p +  labs(x="Maize Area", y="Conflict (Y/N)") +
  geom_smooth(method = "lm") +
  theme(axis.ticks.y=element_blank(), axis.text.y=element_blank()) +
  annotate("text", x = 2, y = 0.25, label = "Correlation = 0.142") +
  geom_point()
dev.off()
###Plot nighttime light
pdf("burntl.pdf")
p <- ggplot(bur.dat, aes(x = nlights_calib_mean, y = civconf))
p +  labs(x="Nighttime Light", y="Conflict (Y/N)") +
  geom_smooth(method = "lm") +
  theme(axis.ticks.y=element_blank(), axis.text.y=element_blank()) +
  annotate("text", x = 0.05, y = 0.25, label = "Correlation = 0.007") +
  geom_point()
dev.off()


###Afghanistan
#Subset data
afg.dat <- subset(state.dat, ccode==700)
##Maps (Figure A10)
#Collapse data
values.a1 <- summaryBy(maize_area ~ latitude + longitude, keep.names = TRUE, FUN=c("mean"), data=afg.dat, na.rm=TRUE)
values.a2 <- summaryBy(nlights_calib_mean ~ latitude + longitude, keep.names = TRUE, FUN=c("mean"), data=afg.dat, na.rm=TRUE)
values.a3 <- summaryBy(civconf ~ latitude + longitude, keep.names = TRUE, FUN=c("mean"), data=afg.dat, na.rm=TRUE)
#Color patterns
col.s1 <- rgb(0,1,0, alpha=1*(values.a1$maize_area/max(values.a1$maize_area)), maxColorValue=1)
col.s2 <- rgb(0,0,1, alpha=1*(values.a2$nlights_calib_mean/max(values.a2$nlights_calib_mean)), maxColorValue=1)
col.s3 <- rgb(1,0,0, alpha=1*(values.a3$civconf/max(values.a3$civconf)), maxColorValue=1)
##Plot
pdf("afgmap.pdf")
#Split to three columns
par(mfrow=c(1,3))
#Maize
map("worldHires", "Afghanistan", col="white", fill=TRUE)
points(values.a1$longitude, values.a1$latitude, pch=15, col=col.s1, cex=1.3)
#Nighttime light
map("worldHires", "Afghanistan", col="white", fill=TRUE)
points(values.a2$longitude, values.a2$latitude, pch=15, col=col.s2, cex=1.3)
#Rebellion
map("worldHires", "Afghanistan", col="white", fill=TRUE)
points(values.a3$longitude, values.a3$latitude, pch=15, col=col.s3, cex=1.3)
dev.off()

#Correlations (Figure A11)
cor(afg.dat$civconf, afg.dat$maize_area)
cor(afg.dat$civconf, afg.dat$nlights_calib_mean)
###Plot maize
pdf("afgmaize.pdf")
p <- ggplot(afg.dat, aes(x = maize_area, y = civconf))
p +  labs(x="Maize Area", y="Conflict (Y/N)") +
  geom_smooth(method = "lm") +
  theme(axis.ticks.y=element_blank(), axis.text.y=element_blank()) +
  annotate("text", x = 0.4, y = 0.25, label = "Correlation = 0.181") +
  geom_point()
dev.off()
###Plot nighttime light
pdf("afgntl.pdf")
p <- ggplot(afg.dat, aes(x = nlights_calib_mean, y = civconf))
p +  labs(x="Nighttime Light", y="Conflict (Y/N)") +
  geom_smooth(method = "lm") +
  theme(axis.ticks.y=element_blank(), axis.text.y=element_blank()) +
  annotate("text", x = 0.1, y = 0.25, label = "Correlation = 0.085") +
  geom_point()
dev.off()

###Colombia
#Subset data
col.dat <- subset(state.dat, ccode==100)

##Map (Figure A12)
#Collapse data
values.c1 <- summaryBy(maize_area ~ latitude + longitude, keep.names = TRUE, FUN=c("mean"), data=col.dat, na.rm=TRUE)
values.c2 <- summaryBy(nlights_calib_mean ~ latitude + longitude, keep.names = TRUE, FUN=c("mean"), data=col.dat, na.rm=TRUE)
values.c3 <- summaryBy(civconf ~ latitude + longitude, keep.names = TRUE, FUN=c("mean"), data=col.dat, na.rm=TRUE)
#Color patterns
col.s1 <- rgb(0,1,0, alpha=1*(values.c1$maize_area/max(values.c1$maize_area)), maxColorValue=1)
col.s2 <- rgb(0,0,1, alpha=1*(values.c2$nlights_calib_mean/max(values.c2$nlights_calib_mean)), maxColorValue=1)
col.s3 <- rgb(1,0,0, alpha=1*(values.c3$civconf/max(values.c3$civconf)), maxColorValue=1)
##Plot
pdf("colmap.pdf")
#Split to three columns
par(mfrow=c(1,3))
#Maize
map("worldHires", "Colombia", col="white", fill=TRUE)
points(values.c1$longitude, values.c1$latitude, pch=15, col=col.s1, cex=1.3)
#Nighttime light
map("worldHires", "Colombia", col="white", fill=TRUE)
points(values.c2$longitude, values.c2$latitude, pch=15, col=col.s2, cex=1.3)
#Rebellion
map("worldHires", "Colombia", col="white", fill=TRUE)
points(values.c3$longitude, values.c3$latitude, pch=16, col=col.s3, cex=1.2)
dev.off()

#Correlations (Figure A13)
cor(col.dat$civconf, col.dat$maize_area)
cor(col.dat$civconf, col.dat$nlights_calib_mean)
###Plot maize
pdf("colmaize.pdf")
p <- ggplot(col.dat, aes(x = maize_area, y = civconf))
p +  labs(x="Maize Area", y="Conflict (Y/N)") +
  #scale_y_continuous(limits = c(0, 1.01)) +
  geom_smooth(method = "lm") +
  theme(axis.ticks.y=element_blank(), axis.text.y=element_blank()) +
  annotate("text", x = 1.5, y = 0.25, label = "Correlation = 0.200") +
  geom_point()
dev.off()
###Plot nighttime light
pdf("colntl.pdf")
p <- ggplot(col.dat, aes(x = nlights_calib_mean, y = civconf))
p +  labs(x="Nighttime Light", y="Conflict (Y/N)") +
  geom_smooth(method = "lm") +
  theme(axis.ticks.y=element_blank(), axis.text.y=element_blank()) +
  annotate("text", x = 0.1, y = 0.25, label = "Correlation = 0.120") +
  geom_point()
dev.off()

##################################
###Descriptive Stats (Table A2)###
##################################

#Rebellion
min(fs.dat$civconf, na.rm=TRUE)
median(fs.dat$civconf, na.rm=TRUE)
mean(fs.dat$civconf, na.rm=TRUE)
max(fs.dat$civconf, na.rm=TRUE)
sd(fs.dat$civconf, na.rm=TRUE)

#Maize area
min(fs.dat$maize_area, na.rm=TRUE)
median(fs.dat$maize_area, na.rm=TRUE)
mean(fs.dat$maize_area, na.rm=TRUE)
max(fs.dat$maize_area, na.rm=TRUE)
sd(fs.dat$maize_area, na.rm=TRUE)

#Nighttime light
min(fs.dat$nlights_calib_mean, na.rm=TRUE)
median(fs.dat$nlights_calib_mean, na.rm=TRUE)
mean(fs.dat$nlights_calib_mean, na.rm=TRUE)
max(fs.dat$nlights_calib_mean, na.rm=TRUE)
sd(fs.dat$nlights_calib_mean, na.rm=TRUE)

#Population
min(fs.dat$pop, na.rm=TRUE)
median(fs.dat$pop, na.rm=TRUE)
mean(fs.dat$pop, na.rm=TRUE)
max(fs.dat$pop, na.rm=TRUE)
sd(fs.dat$pop, na.rm=TRUE)

#Travel time
fs.dat$ttime <- expm1(fs.dat$logttime)
min(fs.dat$ttime, na.rm=TRUE)
median(fs.dat$ttime, na.rm=TRUE)
mean(fs.dat$ttime, na.rm=TRUE)
max(fs.dat$ttime, na.rm=TRUE)
sd(fs.dat$ttime, na.rm=TRUE)

#Rebellion t-1
min(fs.dat$lagcivconflagtemp, na.rm=TRUE)
median(fs.dat$lagcivconflagtemp, na.rm=TRUE)
mean(fs.dat$lagcivconflagtemp, na.rm=TRUE)
max(fs.dat$lagcivconflagtemp, na.rm=TRUE)
sd(fs.dat$lagcivconflagtemp, na.rm=TRUE)

#Rebellion spl.
min(fs.dat$splcivconf, na.rm=TRUE)
median(fs.dat$splcivconf, na.rm=TRUE)
mean(fs.dat$splcivconf, na.rm=TRUE)
max(fs.dat$splcivconf, na.rm=TRUE)
sd(fs.dat$splcivconf, na.rm=TRUE)

#Temperature
min(fs.dat$temp, na.rm=TRUE)
median(fs.dat$temp, na.rm=TRUE)
mean(fs.dat$temp, na.rm=TRUE)
max(fs.dat$temp, na.rm=TRUE)
sd(fs.dat$temp, na.rm=TRUE)

#Drought
min(fs.dat$spi6, na.rm=TRUE)
median(fs.dat$spi6, na.rm=TRUE)
mean(fs.dat$spi6, na.rm=TRUE)
max(fs.dat$spi6, na.rm=TRUE)
sd(fs.dat$spi6, na.rm=TRUE)

#Drugs
min(fs.dat$drug_y, na.rm=TRUE)
median(fs.dat$spi6, na.rm=TRUE)
mean(fs.dat$spi6, na.rm=TRUE)
max(fs.dat$spi6, na.rm=TRUE)
sd(fs.dat$spi6, na.rm=TRUE)

#Polity
min(fs.dat$polity2, na.rm=TRUE)
median(fs.dat$polity2, na.rm=TRUE)
mean(fs.dat$polity2, na.rm=TRUE)
max(fs.dat$polity2, na.rm=TRUE)
sd(fs.dat$polity2, na.rm=TRUE)

#New State
min(fs.dat$nwstate, na.rm=TRUE)
median(fs.dat$nwstate, na.rm=TRUE)
mean(fs.dat$nwstate, na.rm=TRUE)
max(fs.dat$nwstate, na.rm=TRUE)
sd(fs.dat$nwstate, na.rm=TRUE)

#Mountains
min(fs.dat$mnt1, na.rm=TRUE)
median(fs.dat$mnt1, na.rm=TRUE)
mean(fs.dat$mnt1, na.rm=TRUE)
max(fs.dat$mnt1, na.rm=TRUE)
sd(fs.dat$mnt1, na.rm=TRUE)

#Border distance
fs.dat$bdist1 <- expm1(fs.dat$logbdist1)
min(fs.dat$bdist1, na.rm=TRUE)
median(fs.dat$bdist1, na.rm=TRUE)
mean(fs.dat$bdist1, na.rm=TRUE)
max(fs.dat$bdist1, na.rm=TRUE)
sd(fs.dat$bdist1, na.rm=TRUE)

#Cell area
fs.dat$cellarea <- expm1(fs.dat$logcellarea)
min(fs.dat$cellarea, na.rm=TRUE)
median(fs.dat$cellarea, na.rm=TRUE)
mean(fs.dat$cellarea, na.rm=TRUE)
max(fs.dat$cellarea, na.rm=TRUE)
sd(fs.dat$cellarea, na.rm=TRUE)

#Instability
min(fs.dat$instability, na.rm=TRUE)
median(fs.dat$instability, na.rm=TRUE)
mean(fs.dat$instability, na.rm=TRUE)
max(fs.dat$instability, na.rm=TRUE)
sd(fs.dat$instability, na.rm=TRUE)

#Ethnic fractionalization
min(fs.dat$ethfrac, na.rm=TRUE)
median(fs.dat$ethfrac, na.rm=TRUE)
mean(fs.dat$ethfrac, na.rm=TRUE)
max(fs.dat$ethfrac, na.rm=TRUE)
sd(fs.dat$ethfrac, na.rm=TRUE)

#Religious fractionalization
min(fs.dat$relfrac, na.rm=TRUE)
median(fs.dat$relfrac, na.rm=TRUE)
mean(fs.dat$relfrac, na.rm=TRUE)
max(fs.dat$relfrac, na.rm=TRUE)
sd(fs.dat$relfrac, na.rm=TRUE)

#Gems
min(fs.dat$gem_y, na.rm=TRUE)
median(fs.dat$gem_y, na.rm=TRUE)
mean(fs.dat$gem_y, na.rm=TRUE)
max(fs.dat$gem_y, na.rm=TRUE)
sd(fs.dat$gem_y, na.rm=TRUE)

#Oil production
fs.dat$ross_oil_prod <- expm1(fs.dat$logross_oil_prod)
min(fs.dat$ross_oil_prod, na.rm=TRUE)
median(fs.dat$ross_oil_prod, na.rm=TRUE)
mean(fs.dat$ross_oil_prod, na.rm=TRUE)
max(fs.dat$ross_oil_prod, na.rm=TRUE)
sd(fs.dat$ross_oil_prod, na.rm=TRUE)

#Gas production
fs.dat$ross_gas_prod <- expm1(fs.dat$logross_gas_prod)
min(fs.dat$ross_gas_prod, na.rm=TRUE)
median(fs.dat$ross_gas_prod, na.rm=TRUE)
mean(fs.dat$ross_gas_prod, na.rm=TRUE)
max(fs.dat$ross_gas_prod, na.rm=TRUE)
sd(fs.dat$ross_gas_prod, na.rm=TRUE)

